function [fwd_mesh,pj_error] = reconstruct(fwd_mesh,...
						recon_basis,...
                        frequency,...
						data_fn,...
						iteration,...
						lambda,...
						output_fn,...
						filter_n,...
                        wv_array)

% [fwd_mesh,pj_error] = reconstruct(fwd_mesh, recon_basis, frequency,
%                           data_fn, iteration, 
%                               lambda, output_fn, filter_n, wv_array)
%
% Image reconstruction, calls the appropriate reconstruct function
% Do not use for recon with priors or spectral fluorescence.
%
% fwd_mesh is the mesh
% recon_basis is the reconstruction basis (either pixel basis or mesh)
% frequency is the modulation frequency (MHz)
% data_fn is the data, either a variable or filename
% iteration is the max number of iterations
% lambda is the initial regularization value (usually between 1 and 1000)
% output_fn is the output root filename
% filter_n is the number of mean filters per iteration
% wv_array is optional - the array of wavelengths (for spectral mesh)





% If not a workspace variable, load mesh
if ischar(fwd_mesh)== 1
  fwd_mesh = load_mesh(fwd_mesh);
end

% Select femdata program based on mesh.type
if strcmp(fwd_mesh.type,'stnd') == 1        % Use stnd NIRFAST, non-spectral mesh
    if frequency == 0
        [fwd_mesh,pj_error] = reconstruct_stnd_cw(fwd_mesh,...
						recon_basis,...
						data_fn,...
						iteration,...
						lambda,...
						output_fn,...
						filter_n);
    else
  [fwd_mesh,pj_error] = reconstruct_stnd(fwd_mesh,...
						recon_basis,...
                        frequency,...
						data_fn,...
						iteration,...
						lambda,...
						output_fn,...
						filter_n);
    end
elseif strcmp(fwd_mesh.type,'fluor') == 1   % Use stnd fluorfast, non-spectral mesh
  [fwd_mesh,pj_error] = reconstruct_fl(fwd_mesh,...
                                              recon_basis,...
                                              frequency,...
                                              data_fn,...
                                              iteration,...
                                              lambda,...
                                              output_fn,...
                                              filter_n);
elseif strcmp(fwd_mesh.type,'spec') == 1    % Use spectral mesh
    if (nargin < 9)
        wv_array = fwd_mesh.wv;
    end
    if frequency == 0
        [fwd_mesh,pj_error] = reconstruct_spectral_cw(fwd_mesh,...
                              recon_basis,...
                              data_fn,...
                              iteration,...
                              lambda,...
                              output_fn,...
                              filter_n,...
                              wv_array);
    else
      [fwd_mesh,pj_error] = reconstruct_spectral(fwd_mesh,...
                              recon_basis,...
                              frequency,...
                              data_fn,...
                              iteration,...
                              lambda,...
                              output_fn,...
                              filter_n,...
                              wv_array);
    end
end
